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Abstract. The Composite Operator Method (COM) is formulated, its internals 
illustrated in detail and some of its most successful applications reported. COM en- 
dorses the emergence, in strongly correlated systems (SCS), of composite operators, 
optimally deals with their unusual features and implements algebra constraints, and 
other relevant symmetries, in order to properly compute the unconventional prop- 
erties of SCS. 



1.1 Strong Correlations and Composite Operators 

In the last decades, a large part of the research activity in solid state and 
condensed matter physics has been devoted to the study of electronic systems 
with very unconventional properties (i.e. with properties not consistent with 
the Fermi liquid theory). It is commonly believed that the origin of such 
anomalous behaviors should be traced back to the presence of strong electronic 
correlations in such systems. 

Weakly correlated systems are described by Hamiltonians where the 
multi-body terms, and the interactions they describe, feature so little coupling 
constants with respect to the average energies of the one-body terms (e.g. 
wide-band s and p metals) that we can successfully resort to perturbation 
theories. In particular, in such cases, mean-field-like approximations allow to 
diagonalize the Hamiltonian in terms of new independent particles, described 
by canonical operators, that is, operators bounded to satisfy canonical 
commutation relations. On the contrary, in strongly correlated systems 
(e.g. narrow-band d and / metals) [1,2], the interactions are sufficiently intense 
to make completely useless any perturbation scheme that will be just bound to 
fail. In some relevant cases, even more than the strength of the interactions, it 
is their true nature, that is, their actual operatorial form, to bring all troubles 
(e.g. Kondo model [3] ). 



2 



Adolfo Avella and Ferdinando Mancini 



In order to tackle this problem, we need to change perspective and relax 
some of the constraints implied by the quest for new independent particles. 
In particular, we should not expect to be able to correctly describe strongly 
correlated systems only by means of canonical operators. Strong interactions 
modify dramatically the properties of the original particles. What are observed 
are new particles with new peculiar properties entirely determined by the dy- 
namics and by the boundary conditions (i.e. all elements characterizing the 
physical situation under study). These new objects appear as the final result 
of the modifications imposed by the interactions on the original particles and 
embed, since the very beginning, the eflFects of correlations. Thus, the choice of 
new fundamental operators, whose properties will be self-consistcntly deter- 
mined by dynamics, symmetries and boundary conditions, becomes relevant. 

Let us consider a very simple, but very pedagogical, system in order to 
concretely illustrate this occurrence. A lattice of ions sited sufficiently far 
apart to avoid that the electrons can hop from one to another. For the sake 
of simplicity, we can also imagine that only one s-like orbital is active per ion 
and that the Coulomb repulsion is sufficiently intense only between same-site, 
that is, same-orbital, electrons: 

H = -f,Y,cUi)ca{i) + UY, n^^)rn(i) , (1.1) 

i(T i 

where is the chemical potential, i is a vector of the lattice, Ca{i) is the 
destruction operator in the Heisenberg picture {i = {i,t)) of an electron of 
spin (J in a Wannier state centered at the site i, U is the strength of the 
on-site Coulomb repulsion and ncr{i) — cl{i)ccr{i) is the electronic number 
operator. 

Although the Hamiltonian (1.1) looks very simple, there is neither known 
canonical transformation capable to diagonalize this Hamiltonian nor pertur- 
bation scheme, at any order, getting any closer to the exact solution or even 
just grasping its main features. This is simply due to the impossibility to de- 
scribe this system in terms of independent particles represented by canonical 
operators. 

If we analyze the equation of motion of the electronic operator Co-(i) 
d 

i^c„(«) = [ca{i),H] = -nca{i) + Una{i)ca{i) , (1.2) 

we encounter a fermionic operator, that is, an operator constituted by an 
odd number of electronic operators, describing an electron of spin a dressed 
by the charge fluctuations of the electrons of opposite spin (ct). Let us name 
this new operator ridi] = ■Ha(i)ca{i). The equation of motion of r] reads as 

i|^a(i) = [riaii),H] = {U-i,)va(i) . (1.3) 

It is immediate to sec that there exists a companion of 7], which we named 
^aii) = Co-(i) — T]a{i), whosc equation of motion reads as 
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(1.4) 



Equations (1.3) and (1.4) resemble that of independent particles as the cur- 
rent (i.e. the r.h.s of the equation of motion) is directly proportional to the 
operator itself. We define as eigenoperator of a given Hamiltonian, an op- 
erator whose current is just proportional to the operator itself. We will call 
eigenenergy the proportionality constant between the eigenoperator and its 
current (e.g. (JJ — fj) for 77). 

^, together with 77, constitutes a fermionic closed basis for this sys- 
tem. In our definition, a closed basis is the smallest set of operators nec- 
essary to build up the original operators as a linear combination (e.g. c = 
£,+ rf) and which closes the equations of motion. That is, a set of operators 
(Oi, O2) • • • , On) constitutes a closed basis if i§iOp = CpqOq. 

Often, by means of a fermionic closed basis, it is possible to diagonalize 
the related Hamiltonian, that is, to rewrite the Hamiltonian only in terms of 
number-like operators of elements of the operatorial basis. For instance, the 
Hamiltonian (1.1) can be rewritten as 



This diagonalization, although exact, is formal as ^ and 77 are not canonical 
operators and they close commutation relations with their number-like op- 
erators and 77^77 very different from that closed by the electronic operator 
c with n = c^c (see Table 1.1; site indexes have been neglected: they would 
simply lead to Kronecker deltas). 

Equation (1.5) shows that, owing to the presence of the interactions, the 
original electrons c are no more observables and new stable elementary excita- 
tions, described by the operators ^ and 77, appear. The electronic bare energy 
level E = —n splits into two levels (Ei = — /i and E2 = U — fj.) and the 
original electrons turn out to be exactly the worst place to start: no pertur- 
bative scheme leads to the level splitting, which is the main and only feature 
of this very simple system. On the basis of this evidence, one is induced to 
move the attention from the original electronic operator to the new operators 
generated by the interaction. Such new operators can be written in terms of 
the original ones and are known as composite operators. A formulation 
which would treat composite operators as fundamental objects, instead of the 




(1.5) 
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original electronic operators, would be very promising. As a matter of fact, 
such a formulation would allow to set up approximation schemes where some 
amount of the interaction is already exactly taken into account through the 
chosen operatorial basis, and would permit to overcome the hoary problem of 
finding an appropriate expansion parameter. 

However, one price must be paid. In general, composite operators are nei- 
ther Fermi nor Bose operators, since they do not satisfy canonical (anti)commutation 
relations (see, for instance, Table 1.1). This very simple fact makes a tremen- 
dous difference with respect to the case in which one deals with the original 
electronic operators, which satisfy a canonical algebra. In developing per- 
turbation schemes where the building blocks are propagators of composite 
operators, one cannot use any more the consolidated scheme: diagrammatic 
expansions, Wick's theorem and many other standard tools are no more valid 
or applicable. The formulation of the whole Green's function method must 
be revisited [4, 5] and new frameworks of calculations have to be formulated. 
In the following section, we will formulate the Composite Operator Method 
(COM) [5] and illustrate its internals; COM systematically endorses the emer- 
gence, in strongly correlated systems, of composite operators and optimally 
deals with their unusual features. 



1.2 The Composite Operator Method (COM) 
1.2.1 Basis 

Given a Hamiltonian H describing a solid state system, first we have to 
choose a basis ip constituted of a finite number of composite operators, either 
fermionic or bosonic. A fermionic basis will be used if the original interacting 
particles are electrons and we want to analyse thermodynamic and single- 
particle properties of the system. A bosonic basis can be either the natural 
basis for systems constituted by bosons (e.g. spins, phonons, . . . ) or a set 
of bosonic operators, constituted by an even number of electronic operators, 
representing the electronic number operator (charge), the electronic spin op- 
erator, . . . and necessary to study the charge, spin, . . . response of the system. 
In this latter case, the dynamics described by the bosonic basis is strongly 
coupled to that described by the related fermionic basis and viceversa. As a 
matter of fact, this is the case even if we do not explicitly open the bosonic 
sector as the self-consistent parameters appearing in the fermionic sector are 
just thermal averages of the related bosonic operators. A fully simultaneous 
self-consistent solution of both sectors is usually required. This consequence of 
the non-canonical algebra satisfied by composite operators is a manifestation 
of the relevance (actually, the dominance) of spin, charge, . . . correlations in 
all properties exhibited by strongly correlated systems. 

Obviously, a closed basis is the best basis one can adopt. It is possible 
to obtain a closed basis just crossing the entire hierarchy of the equations 
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of motion of the relevant canonical operator (e.g. the electronic operator, its 
number or spin operator, ...). Namely, one has to: (i) compute high-order 
time derivatives of the relevant canonical operator by repeatedly commuting 
it with the Hamiltonian; (ii) acquire in the basis all distinct composite oper- 
ators appearing in the high-order currents; (iii) terminate as soon as no new 
composite operator appears. This well-defined procedure, being equivalent to 
compute all spectral moments (i.e. the moments of the relevant canonical oper- 
ator Green's function) [6 8] , catches, in principle, all scales of energy featured 
by the system under analysis. 

Following this procedure, one can see that there is a large class of systems 
(finite systems [4,9], bulk systems with interacting localized electrons [10 12], 
Ising-like systems [13-16], . . . ) where it is possible to obtain a closed basis with 
a finite number of elements and one can aim at the exact solution. In all other 
bulk systems, this procedure will lead to a closed basis with an infinite number 
of components. In such cases, one can: rewrite H as H = Ho + Hj, where Hq 
is the most relevant part whose finite closed basis is known; adopt the closed 
basis of Hq as a truncated basis for H; treat Hj as a perturbation to Hq. If 
the system under study features only analytical scales of energy (i.e. scales 
of energy whose mathematical expressions are expandable in power series of 
Hamiltonian model parameters), or if the range of temperatures/frequencies 
of interest does not involve not-analytical scales of energy, the adoption of 
a truncated basis can lead to an accurate description of the system under 
analysis. Such a description can be systematically improved adding to the 
truncated basis more and more elements according to the above described 
procedure. For instance, in the Hubbard model, we need at least a 2-element 
basis to correctly describe the Hubbard on-site Coulomb repulsion U and 
at least a 4-element basis to correctly describe the virtual exchange energy 
J = U'^/U [9,17]. 

Not analytical scales of energy can be present in some models. For example, 
in the single-impurity Anderson model [3], the emergence of a not-analytical 
Kondo energy scale is directly connected to the virtual exchange processes 
triggered by the hybridization between conduction and impurity electrons and 
dominated by on-site Coulomb repulsion at the impurity site. This occurrence 
forces to explore different routes to construct a proper operatorial basis. In 
fact, in such cases, the adoption of any tnmcated basis would simply lead to 
a description of the system under analysis completely lacking any reference to 
not-analytical scales of energy. Such a description would turn out really poor 
at low enough temperatures/frcqiiencies where neglected energy scales usually 
manifest themselves triggering quite unconventional behaviours (impurity spin 
screening, superconductivity, . . . ). Actually, in such cases, one has to properly 
tweak the equations of motion of specific composite operators and smartly 
exploit self-consistency in order to get sensible results by means of a basis 
with a finite number of elements [18, 19]. 

Summarizing, many recipes can assure a correct and controlled description 
of relevant aspects of the dynamics. One can put in the basis: 
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• the higher-order composite operators emerging from the hierarchy of the 

equations of motion: the conservation of a definite number of spectral 
moments is guaranteed; polynomial and analytical scales of energy can be 
resolved with desired accuracy [5,8] 

• the cigcnopcrators of Hamiltonian terms describing specific interactions or 
virtual processes: the latter are correctly taken into account [9, 20] 

• the eigenoperators of the problem reduced to a small cluster: all polynomial 
and analytical energy scales and their competitions/interplays are just 
exactly taken into account at short distances/large momenta [9,21-23] 

• the composite operator describing the Kondo-like singlet [18, 19, 24, 25] 

1.2.2 Equations of motion of ip 

Once the basis ip{i), constituted of n composite operators, has been chosen, 
we adopt a matricial notation and write 



/Mi)' 

m = ] 



(1.6) 



where we have not specified the nature, fcrmionic or bosonic, of the basis. 
In case of fermionic operators, it is understood that we use the spinorial 
representation 

f^-" 

The current J of tp can be exactly rewritten as 

= i^ v(i) = mi), H] = Y1 ^(i' j)V'a. t) + , (i.s) 

j 

where the matrix e is known as the energy matrix. 6J is named residual 
current and describes the component of the current J orthogonal to the basis 

[6J{i,t),i;Hi,t)]^) = ^ e(i,j) = ^m(i, 1)7-1(1, j) . (1.9) 

1 

It is absolutely worth noticing that iff is a closed basis, SJ is identically 
zero. 

m(i,j) = ^[J(i,t), V'^(j,i)]^^ is simply named m-matrix and /(i,j) = 
[tp{i, t), "(/"^(j, t)] ) is known as the normalization matrix. The matrix / is 



V/ 

hermitian by definition, while the matrix m is hermitian because of the time 
independence of the matrix /: = i^7 = m—rri^. {■■■) stands for the quantum 
mechanical average in the (grand-)canonical ensemble. 77 = 1, —1 marks anti- 
commutators {,} and commutators [, ], respectively. For fermionic (bosonic) 
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composite operators, the choice 77 = 1 (77 = —1) is much more convenient (i.e. 
it makes the calculations much simpler), although the opposite choice is also 
possible. 

Once a basis has been chosen, either closed or truncated, we need to calcu- 
late the normalization / and the m matrices in order to construct the energy 
matrix e, which will soon turn out to be a fundamental quantity. The cal- 
culation of the relevant (anti)commutators is straightforward, but it can be 
lengthy and ciimbcrsome if the basis contain many elements. The thermal av- 
eraging procedure is instead much less mechanical as it involves choosing the 
various phases to be studied. All correlation functions present in the normal- 
ization I and in the energy e matrices need to be computed self-consistently 
and form a first block of unknowns in the theory. 

Weight and orthogonality of composite operators 

The entries of the normalization matrix / give a measure of both the weight 
of a composite operator (diagonal entries) and the degree of orthogonality 
between two of them (off-diagonal entries). These two concepts (weight and 
orthogonality) have obvious meanings for canonical electronic operators whose 
anti-commutation relations are just C- numbers ({c;,C;,} = Sw): an electronic 
operator always weights just 1 and is orthogonal by definition to any other 
electronic operator, clearly describing a single-particle state with different 
quantum numbers. As for composite operators, owing to the non-canonical 
commutation relations they obey, these concepts become more relevant and 
less obvious: both properties strongly depend on the actual operatorial form 
of the composite operators as well as, through the thermal averaging, on the 
Hamiltonian describing the system under analysis and the boundary condi- 
tions. Such a dependence manifests itself through the emergence, in the / 
matrix entries, of thermal averages of bosonic composite operators, whose de- 
termination is not always straightforward. Actually, it is just the presence of 
these parameters in the / matrix, to leave sometimes a problem unsolved al- 
though a closed basis has been individuated. If the basis is not closed {SJ ^ 0), 
the same kind of averages appears in the e matrix too. 

These features of composite operators (non-trivial weights and imperfect 
orthogonality) are distinctive landmarks of strong correlations. They are fun- 
damental to the occurrence of well-known effects such as the transfer of spec- 
tral weight (directly connected to composite operator weights) between (sub- 
)bands on changing model or external (temperature, filling, pressure, external 
fields, . . . ) parameters, and the interplay among scales of energies (connected 
to the overlapping of composite operators) . 

For instance, if we choose ip = (^, 77) as (closed) fermionic basis for the toy 
model (1.1), we would find the following I and e matrices: 
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'(iJ) 



1 - (nff(i)) 
(n^(i)) 

-H 



1 
-1 



SiAa' I ?7 - /X 



(1.10) 



As the chosen basis is closed, the energy matrix is free of unknown param- 
eters, while the normalization matrix contains anyway some of them to be 
computed self-consistently according to the chosen phases to be investigated 
(e.g. paramagnetic, ferromagnetic, Neel, CDW, SDW, ...). 



1.2.3 Dyson equation for G 

We can now define a generalized matricial Green's function G for the basis 



G^{i,j) = (Q[m^Hj)W 



e{u-t,){m),^^{j)],) for Q = n 

~e{t,~U){[i;{z),^^Hj)]v)ior Q = A 

(1.11) 

where TZ, A and C stand for the usual retarded, advanced and causal operators, 
respectively. 

For the sake of simplicity, we pick up a phase with spatial homogeneity and 
move to momentum and frequency spaces through Fourier transform. Then, 
according to (1.8), G satisfies the following equation of motion 

G2(k,a;) = G2(k,a;) + GC(k,a;)T2(k,a;)G2(k,a;) , (1.12) 

where Gq (k, uj) , determined by the equation 

[w-£(k)]Go2(k,c.)=/(k) , (1.13) 

will play the role of fundamental building block and the scattering matrix 
T^(k, w) is defined as 

T^iKcu) = /-i(k)B2(k,w)/-i(k) , (1.14) 

withBQ{i,j) = {Q[6J{t)6J^{j)]^). 

Now, we can introduce the irreducible self-energy U^{k, w), defined through 
the relation T2(k,a;)G^(k,a;) = I-\k)E^{k,uj)G^{k,uj), and get 

GQ(k,a;)= ,, \on...^ Ii^) 



w-e{k)-IJQ{k,uj) ' ' [G2(k,a;)]-i-/-i(k)i:2(k,w) ■ 

(1.15) 

Equation (1.15) generalizes the well-known Dyson equation to the case of 

Green's functions of composite operators [4,5]. Together with the apparent 
similarities, it is absolutely worth noting the striking differences: 
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• all ingredients are matrices either n x n (bosonic) or 2n x 2n (fermionic); 

• the equivalent of the usual non-interacting Green's function, Gq, is, in- 
stead, fully interacting and, in principle, exact (if the basis tp is closed 
then we have 5 J = 0, consequently S = 0, and finally G = Go); 

• the normalization matrix / evidently assumes the role of overall weight of 
the matricial Green's function G (/ is its 0*'' moment: lim^^-^oo G ~-> -f/w); 

• the eigenvalues of the energy matrix s clearly play the fundamental role of 
poles, which will be more or less severely modified by the self-energy E. 

• the self-energy S describes both the dynamics orthogonal to that described 
by the chosen basis tp, and the residual interactions among the elements 
of the latter. Smarter is the choice of tjj, up to a closed basis, less relevant 
is the knowledge of S to fully understand the system under study. Hence, 
we decided to rename the generalized irreducible self-energy S as residual 
self-energy: it is just the irreducible propagator of the residual current 6 J. 

1.2.4 Propagator Go 

Let us take a step back, the propagator Go satisfies the following equation 

[oj-e{k)]G^ik,w)=I{k) , (1.16) 
whose most general solution is [4,5] 

Gf{Kuj) = Y,\r 
1=1 ^ 

(T^')(k) are the spectral density functions in matrix form, fully determined by 
the energy matrix e and the normalization matrix / as 

(k) = E ^/c' (k)/c6(k) , (1.18) 

c 

where /2(k) is the matrix whose columns are the eigenvectors of the energy 
matrix e. u}i{'k) are the eigenvalues of this latter, g^^^^ik) are unknown mo- 
mentum functions, in matrix form, not fixed by the equations of motion (i.e. 
they correspond to constants in time), to be determined taking into account 
the boundary conditions specific of the type of propagator under study (re- 
tarded, advanced or causal). 

For both fermionic (t? = 1) and bosonic {r] = —1) basis, the boundary 
conditions turn out to be sufficient to fully determine ^'('^^'•^(k) and we obtain 
a Lindhard-like representation for G^"^{'k,Lj) 

Contrarily, ^('^''(k) and, consequently, the causal Green's function Go(k, lo) 
can be fully determined by the boundary conditions only for fermionic basis 



a^'J(k) 
W - W;(k) 



i7r5[a;-a;i(k)]5«S(k)| . (1.17) 
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Gg(k,u;) = Va«(k) 



1=1 



1-/fH ^ /fH 



LO — LOl{k.) + 16 L0 — LJl{k.) — i6 



(1.20) 



where /f('^) is the Fermi distribution function. 

We conclude the discussion about fcrmionic basis reporting the expression 
for the correlation function C{i,j) = (V'(*)V'^(j)) 

n 

C(k,a;) = 2^^[l-/F(a;j(k))]aW(k)<5[c.-a;,(k)] . (1.21) 
1=1 

It is necessary noticing that the correlation function C can be computed by 
means of (1.21) only if the basis tp is closed or if we neglect the residual 
self-energy. Otherwise, we should use the more general expression 

C(k,w) =T2[l-/F(c<;)]9[G^-^(k,a;)] . (1.22) 

For bosonic basis, instead, g*^''''(k) cannot be fully determined by the 
boundary conditions, but we can still obtain a modified Lindhard-like rep- 
resentation for Go(k, w) [4,5] 



Gg(k,w) = -27rir(k)5(w) + ^ (T^'Hk) 



1=1 



l + fB{to) /b(w) 



_u) — u) I {]<.)+ iS u) — u)i{]s.) — iS _ 

(1.23) 

where /b(w) is the Bose distribution function, provided that we explicitly 
keep in the final expression an unknown zero frequency momentum function 
in matrix form -r(k). 

The actual value of r{k) is directly related to the degree of ergodicity of 
the dynamics of the basis ip driven by the Hamiltonian H 

r(i - j) = lim i r dt (V(i, 0)V^ {It)) . (1.24) 

If we would know, by any source of information, that the dynamics of ij) driven 
by H is fully ergodic, the last equation would simply give 

r(i-j) = (V'(i))(Vt(j)) , (1.25) 

leading to an effective removal of the unknown F function from the theory. 

Unfortunately, although many people just overlook this deliberately, there 
is absolutely no way to asses ergodicity in advance and -r(k) forms another 
block of unknowns in the theory. In particular, ergodicity cannot be supposed 
a priori for finite systems treated in statistical cnscimbles different from the 
micro-canonical one and for bosonic operators commuting with the Hamilto- 
nian (i.e. for constants of motion). 

r{k) also appears in the expression for the correlation function C{i,j) = 
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(7(k,w) =27rr(k)5H + 27r^[l + /B(wz(k))]aW(k)5[a;-a;,(k)] , (1.26) 

1=1 

but this is absolutely not the only issue with the above expression (1.26). 

If C7(')(k) is not identically zero for all values of k where Wi(k) vanishes, 
the correlation function C(k, w) diverges as l//3ti;((k) for the same values of 
k in the limit w — >■ (i.e. at all times). Such a behaviour of C(k, w) usually 
manifests the establishing of long-range spatial correlations in the system, but 
it is admissible iff the corresponding correlation function in real space C(r, t) 
stays always finite. Accordingly, the divergence should be integrable; this im- 
mediately excludes the possibility to have long-range order in finite systems 
and, at finite temperatures, in infinite systems (i.e. in the thermodynamic 
limit) with too low spatial dimension. Actually, the lowest spatial dimension 
allowed to host long-range order will simply depend on the actual functional 
form of the vanishing (jJii}i) (e.g. if Wi(k) oc |k|" then the spatial dimension has 
to be strictly larger than a). No restriction applies to infinite systems of any 
dimension at zero temperature. This is just the content of Mermin- Wagner 
theorem [26]. 

For bosonic basis too, it is necessary noticing that the correlation function 
C can be computed by means of (1-26) only if the basis ip is closed or if we 
neglect the residual self-energy. Otherwise, we should use the more general 
expression 

C{\.M = -'^ lllff, ^ [G'iKo^)] . (1.27) 

J- + '^JB\^) 

The use of the casual Green's function, for a bosonic basis, is strictly neces- 
sary in order to properly take into account the ergodicity issue that does not 
manifest in neither the advanced nor the retarded propagators, but only in 
the causal propagator and in the correlation function. 

1.2.5 Residual self-energy S 

Obviously, we do not need to compute the residual self-energy at all (it is just 
identically zero) if we choose a closed basis; this is one of the main reasons 
why a closed basis is the best one you can choose. Accordingly, as much as 
a truncated basis is large enough less relevant will be the contribution of the 
residual self-energy to the description of the system under analysis. It is worth 
noticing that even if we would completely neglect the residual self-energy, the 
Green's function of the original interacting particles constituting the system, 
expressed in terms of the relevant entries of the propagator Go, will anyway 
feature a fully momentum and frequency dependent irreducible self-energy 
with a (n — l)-polar structure, but sufficient to describe all scales of energy 
caught by the chosen basis 

Nevertheless, one cannot neglect the residual self-energy without being 
aware of the main drawback: one is actually promoting to the rank of true 
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particles (i.e. with an infinite life-time) objects that, being still subject to 
all virtual processes not properly taken into account by a truncated basis, 
have finite life-times (i.e. they are quasi-particle) roughly inversely propor- 
tional to the largest neglected energy scale involving them (i.e. the transition 
described by a composite operator can or cannot be one of those necessary 
to construct the relevant virtual process). Clearly, this can be systematically 
controlled by enlarging the basis, although only on a quantitative level. Again, 
not-analytical energy scales can change so dramatically the properties of par- 
ticles and quasi-particles to require a qualitative change of perspective also 
regarding the residual self-energy determination and role, as described in the 
Sect. 1.2.1. 

As a matter of fact, it is sometimes convenient keeping the operatorial 
basis somewhat simpler than what is actually handleable in order to get sim- 
pler expressions for the residual currents too and effectively compute, starting 
from the latter, the residual self-energy. Such a procedure can lead to a de- 
scription of the system under analysis featuring both some of the relevant 
energy scales and the decay effects inherent to the quasi-particle nature of 
composite operators belonging to a truncated basis. 

In these years, we have been developing and testing different ways to 
compute residual self-energy; among others: the Two-site resolvent approach 
[27, 28] and the Non Crossing Approximation [29-36] . Any of them has its 
pros and cons and specific problems can be better tackled by one or the other 
(see Sect. 1.3). 

1.2.6 Self-consistency 

The intrinsic complexity of the operatorial algebra closed by composite op- 
erators (for instance, see Table 1.1) hides a noteworthy possible exploitation 
of the same algebra. Composite operators, whose lattice sites overlap (any 
composite operator may span a certain number of lattice sites), obey algebra 
constraints directly coming from the Pauli exclusion principle. For instance, 
^ and T] satisfy the following exact relation: ^cr{i)'nt'i^) ~ ^ that can be easily 
traced back to Co-(i)ccr(i) = 0. 

The relevance of such algebra constraints resides, on one side, in their 
capability to enforce the Pauli principle and its derivatives, and, on the other 
side, in the not-so-trivial request that they should be obeyed at the level of 
thermal averages too. In fact, the related thermal averages (e.g. (^CT(i)??J./(i)) = 
0) , through the well-known relation existing between correlation functions and 
Green's functions (fluctuaction-dissipation theorem; see Sect. 1.2.4), depend 
on all unknown parameters appearing in the relevant Green's function (in 7, 
e, and S plus r'(k) for bosonic basis) and, in turn, can be used to fix them: 

(V(i)V^(j)) = C(i-j,t = 0). (1.28) 

The l.h.s of (1.28), for i and j such that elements of the basis ip span on, 
at least, one common site, would be fixed by the algebra, which imposes 
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contractions (e.g. Ccr(i)^I'(i) = Oj n(i)n(i) = n(i) + 277|(i)?7-f (i)). The r.h.s. of 
(1.28) is given by the actual expression of the Green's function and contains 
ah unknowns of the theory (for bosonic basis, it also contains r{i — j)). 

This deceptively simple conclusion has enormous implications on both the 
capability to solve, either exact or approximately, strongly correlated systems, 
and the quality of the solution. Not only algebra constraints allow to find a 
solution, but they make this solution as closer as possible to the exact one 
because they embody the primary cause of electronic correlations: the Pauli 
principle. 

It is also worth noticing that the symmetries enjoined by the Hamiltonian 
imply the existence of constants of motion and the possibility to formulate 
relations among matrix elements of the relevant Green's functions known as 
Ward-Takashi identities [37,38]. These identities can/should also be used to 
fix the unknowns and to constrain the theory. 

1.2.7 Summary 

Summarizing, COM framework envisages four main steps: 

1. Choose a composite operator basis according to the system under analy- 
sis and all information we can gather from relevant numerical and exact 
solutions (see Sect. 1.2.1) 

2. Compute / and m matrices and obtain e matrix and propagator Go in 
terms of unknown correlators and unknown F function (see Sect. 1.2.1 
and 1.2.4). 

3. Choose a recipe to compute S or just neglect it (see Sect. 1.2.5). 

4. Self-consistently compute the unknowns through Algebra Constraints and 
Ward-Takashi identities (see Sect. 1.2.6). 

In the last fifteen years, COM has been applied to several models and mate- 
rials: Hubbard [4,5,9,17,28,32,36,39-42], p-d [43-45], t-J [9,27], t-t'-U [46-48], 
extended Hubbard (t-U-V) [49,50], Double-exchange [22], Kondo [18], An- 
derson [19], Kondo-Heisenbcrg [51], two-orbital Hubbard [52], Kondo lat- 
tice [53], Ising [10-14,16], BEG [54], Heisenberg [55], J1-J2 [56-58], Hubbard- 
Kondo [59], Singlet-hole [60], Cuprates [31,61-63] . . . A comparison with the 
results of numerical simulations has been systematically carried on. The in- 
terested reader may refer to the works cited in Ref. [5] and, for the last years, 
at the web page: http://scs.physics.unisa.it. 

In the second part of this manusc;ript, as relevant application of the formal- 
ism, we will consider the Hubbard model and we will go through the different 
approximation schemes illustrated in the previous Sections in a systematic 
way. A comprehensive comparison with the results of numerical simulations 
and with the experimental data for high T^, cuprates will be also reported. 
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1.3 Case study: The Hubbeird model 

1.3.1 The Hamiltonian 

The Hubbard model reads as 

H = J^ (-M^ij - 2dtaij) cHi)c{j) + C/^nt(i)n4(i) . (1.29) 

ij i 

The notation is the same used in Sect. 1.1 with the following additions: t is 
the hopping and the energy unit; d is the dimensionality of the system; ay is 
the projector on the nearest-neighbour sites, whose Fourier transform, for a d- 
dimensional cubic lattice with lattice constant a, is Q;(k) = g X]n=i cos(fc„a). 
The electronic operators c(i) and c^(i), as well as all other fermionic operators, 

are expressed in the spinorial notation [e.g. c^(i) — ^c|(i) c|(z)^ . A detailed 

and comprehensive summary of the properties of Hamiltonian (1.29) are given 
in Sect. 1 of Ref. [5]. We here report a study of this model by means of the 
Composite Operator Method as formulated in Sect. 1.2. In order to proceed 
in a pedagogical way, we will go through different stages. In Sect. 1.3.1, we 
consider a truncated fermionic basis (2-pole approximation) given by the two 
Hubbard operators ^(i) and r](i). In Sect. 1.3.1, we complement the analysis 
considering the related bosonic sector (charge and spin). In Sect. 1.3.2, we im- 
plement the study by considering the contribution of the residual self-energy. 
In Sect. 1.3.3, we enlarge the truncated basis by including higher-order com- 
posite fields (4-pole approximation). At all these stages, we will present COM 
results and compare them with experimental and/or simulation data. Due to 
the pedagogical nature of this manuscript, we will restrict the analysis to the 
paramagnetic state. Ordered phases (ferro and antiferromagnetic phases) have 
been also analyzed and the related results can be found in Ref. [5] . Finally, in 
Sect. 1.3.4, we will discuss the superconducting solution of the model in the d- 
wave channel and compare COM results with high-Tc cuprates experimental 
data. 



Two-pole solution - Fermionic Sector 

On the basis of the discussions reported in Sects. 1.1 and 1.2.1, we will adopt 
as fermonic basis 




(1.30) 



where and ri{i) are the Hubbard operators defined in Sect. 1.1. This field 
satisfies the equation of motion 

dt ^ ^ \{U - fj,)r]{i) + 2dtTr{i) J ' ^ ^ 
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with n^(i) = c'{i)afiC{i) and 7r(i) = ^a'^nfi{i)c"'{i) + c{i)c'°'{i)c{i). = 
(— l,cr) and ct^ = (l,cr), where cr arc the Pauh matrices. Hereafter, given 
a generic operator ^(i), we will use the notation = X^j ciij^O, i)- The 

current J{i) is projected on the basis (1.30) and the residual current SJ{i) is 
neglected (truncated basis). The normalization and energy matrices have the 
expressions 



hi 0\ \ _ /mn(k)/iVTOi2(k)/. 



^22 



0/22; V n/2; ^ ^ V"il2(k)-rri "J22(k)/22 

where the entries of the m-matrix are given by 



1 ™ n,\ T-i 
22 . 

(1.32) 



mii(k) = -/i/ii - 2dt [A + a(k)(l - n +p)] (1.33) 

mi2(k) = 2dt[Z\ + a(k)(p-/22)] (1.34) 
TO22(k) = {U- m)/22 - 2dt [A + a(k)p] , (1.35) 

n = 1/-^ Si is the particle number per site; the parameters A and p 
cause a constant shift of the bands and a bandwidth renormalization, respec- 
tively, and are defined as 

A = {Cm\^))-{v''i^)vH^)) (1-36) 
f =^KWn,.W)-([ctWc4W]"c|(*)4(i)). (1.37) 

A is the difference between upper and lower intra-subband contributions to 
kinetic energy, while p is a combination of the nearest-neighbor charge-charge, 
spin-spin and pair-pair correlation functions. 

The retarded G^{i,j) = ('7^ {V'(^)V'^(i)}) and the correlation C{i,j) = 
(V'(»)V'^(i)) functions are given by 

n=l ^ ' 

2 

C(k,u;) = [1 - /fM] (7(")(k)<5[c. - £^„(k)] . (1.39) 

n=l 

The energy spectra -B„(k) read as 

£;i(k) = i?(k) + Q(k) i;2(k) = i?(k) - Q(k) , (1.40) 

where 



i?(k) = 1 [C/ - 2M - 4.to(k)] - ^ Q(k) = ll^(k) + 1^ , 

11 22 y 11 22 

(1.41) 
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with gCk) = —U + /^^/"^ '^i2(k)- The spectral functions ^("'^(k) have the 
following expressions 

(k) = ^ [1 + ^] (k) = ^ [1 - sgg) ] 

^'Hk) = l^ a['^{k) = -W!§ ■ (1-42) 



'12 y^J — 2Q(k) "12 \'^J — 2Q(k) 

^^22 — 2 l-*- 2Q(k)J "22 K'^J — 2 l"*" ^ 2Q(k) 



The determination of G'^(k, a;) and C(k,w) requires the knowledge of 
the chemical potential and of the two bosonic correlators A and p. These 
quantities are self-consistently determined by means of the following set of 
coupled equations 

n = 2(1 - Cii - C22) (1.43) 
A = Q\ - (1-44) 
C,2 = mvHi))=0, (1.45) 

where Cab = (V'a(i)V'6(i)) a-nd C^;, = (i)'!/'5(i))- The first equation fixes the 
chemical potential in terms of all other parameters; the second comes from the 
definition of A [cfr. (1.37)]; the third comes from the constraint (1.28). Once 
this set of coupled self-consistent equations has been solved, we can calculate 
all relevant single-particle and thermodynamical properties of the model. The 
double occupancy per site: D = 1/N '^i{n-f{i)n^{i)) = I22 — C22] the energy 
bands £'n(k); the Fermi surface by means of the equation ^^^(k) = 0; the 
momentum distribution function n(k) and the density of states N{uj): 



/+00 
-00 



{2n) 



Ob 



TT 

--9[GS(k,a;) 

TT 



(1.46) 
(1.47) 



where 17b is the volume of the first Brillouin zone, G^(k, cj) = G'^(k, u); 

the internal energy E = (H) = 4dt j,^^ C"^ + UD; the specific heat C = 
dE/dT; the free energy F(n, T) = fi{n' , T)dn'; the entropy S = {E-F)/T. 

Once the fermionic propagator is known, there are several ways to compute 
response functions (i.e. the retarded propagators of the two-particle excita- 
tions: charge, spin, pair, . . . ). These techniques are usually based on diagram- 
matic expansions of the two-particle propagators in terms of the single-particle 
one. However, when operators with non-canonical commutations are involved, 
the complicated algebra invalidates the Wick theorem and, consequently, does 
not allow any simple extension of decoupling schemes and more involved dia- 
grammatic approximations [64,65] are needed. Another technique [5,66], the 
one-loop approximation for composite operators, has been developed by means 
of the equations of motion approach. By using this technique it is possible to 



1 The Composite Operator Method (COM) 



17 



calculate the causal function (T[n^(i)n^(j)]) and then the charge and spin 
susceptibilities. We do not report the details of the calculations, which can be 
found in Ref. [5], and summarize the results. For the causal propagators we 
have 



a,6,c— 1 



,,b,c=l 



where Q'^abcA^'-j) = G'^b(''i the fermionic loop constructed from 

the propagator G'^^{i,j) = (T[f/'a(*)V'b(j)])- ^'-'^ charge and spin suscepti- 
bilities, we have 

^^^''''^)='n^ ^ 4-;ebac(k,^) (1.50) 

a,fc,c— 1 

^-('^'"^= n + 2D-l^ ^ CQS.c(k,-), (1.51) 

a,b.c—l 

where Q^^^(k, a;) is the retarded part of (5^^^^(k, a;) and has the expression 
' r .{MEM - /F[i^n(k + P)]}ai:)(k + p)a(™)(p) 



^ n,m=l 



Results and comparisons 



(1.52) 



Algebra constraints The violation of the Algebra Constraints in the Roth 
and Hubbard I approximations is analyzed in Fig. 1.1 (left panel), where 
the normalized Pauli amplitude Ap = C12/C22, which is identically zero in 
COM (as it should be due to the Pauli principle), is reported versus n at 
T = for various values of U. We see that in the Roth [68, 69] and Hubbard 
I [1] schemes the Pauli principle is satisfied only in the cases n = and 
n = 1. The deviation increases by decreasing U, reaching its maximum value 
in the noninteracting limit (i.e. U = 0). On the contrary, the Pauli principle is 
recovered in the limit ?7 — )■ 00 for any value of n. At any rate, the Hubbard I 
solution violates several othcir sum rules [70] . As a consequence of the fact that 
the Pauli principle is not satisfied, the particle-hole symmetry enjoined by the 
Hubbard Hamiltonian is also violated [5,71]. In order to give a measure of 
how relevant is fulfilling Algebra Constraints, and in particular those directly 
coming from Pauli principle, in the right panel of the same Fig. 1.1, we show 
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Fig. 1.1. (left) Pauli amplitude Ap — C12/C22 versus filling n at T = for several 
values of [/. (right) Double occupancAr D versus filling n for T = 1/6 and U = 2. The 
solid, dotted and dashed lines stand for COM, Hubbard I and Roth, respectively. 
qMC data refer to a 12 x 12 cluster [67]. 




-4 -3 -2 -1 1 2 3 4 1 2 3 4 5 

h T 



Fig. 1.2. (left) Filling n versus chemical potential for T = 0.05, 1, 2 and 4 

and U — 8. BA result is taken from Ref. [72]. (right) Chemical potential /i versus 
temperature T for f/ = 4 and n = 0.5 0.95. FTLM result is taken from Ref. [73]. 



COM results for the double occupancy in comparison with the Roth and 
Hubbard I ones, and the data obtained on finite size clusters by quantum 

Monte Carlo method [67]. 

Chemical potential In Fig. 1.2 (left panel) , we show the particle density versus 
chemical potential ^ for various temperatures. COM results are compared 
with the Bethe Ansatz (BA) ones [72]. The agreement is very good at low 
temperatures for densities smaller than 0.55. In the half-filled chain, T « t is 
a relevant temperature as it signs the border between T-regions dominated 
by either spin or charge correlations [74,75]. The agreement between COM 
and BA at T > Hs very good for the whole range of filling. Of course, at high 
temperatures, COM result reaches an excellent agreement with BA since the 
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Fig. 1.3. (left) Double occupancy D versus Coulomb interaction (7 at T = and n = 
1. Dotted lines stand for BA results, (right) Double occupancy D versus Coulomb 
interaction Z7 for T = 1/6 and n = 1, sohd symbols stand for qMC data (4 x 4) [76]. 



effect of correlations is completely suppressed. In Fig. 1.2 (right panel), we 
show COM results for the chemical potential as a function of the temperature 
T and of the filling for J7 = 4 in comparison with numerical results obtained 
by means of the Lanczos technique on a 4 x 4 cluster [73]. The agreement 
is rather good in particular for low filling and high temperature where the 
numerical results are more reliable and the finite size of the cluster is less 
effective. 

Double occupancy The double occupancy for the one-dimensional case is 
shown in Fig. 1.3 (left panel), where COM is plotted versus Coulomb interac- 
tion f/ at T = and for n = 1. The results are compared with the exact ones 
by Bethe ansatz. The agreement with BA is excellent. Such a good agreement 
is not reached by any other analytical approach, neither by the Gutzwiller or 
the ladder and the self-consistent ladder approximations [77-79]. In particu- 
lar, these approaches fail to reproduce the correct asymptotic behavior. As 
shown in Fig. 1.3 (left panel), the double occupancy goes to zero as C/ — >■ oo: 
the electrons localize only at infinite U , where the half-filled Hubbard chain 
is equivalent to the spin- 1/2 AF Heisenberg chain. The double occupancy for 
the two-dimensional case is shown in Fig. 1.3 (right panel), where COM so- 
lution is plotted versus Coulomb interaction 17 at T = and for n = 1. The 
results are compared with the data of numerical simulation by qMC [76] . The 
agreement is excellent. 

Internal Energy The internal energy at n = 1 and T = is shown as a func- 
tion of U in Fig. 1.4 (left panel). The results obtained by means of the BA [80] 
and other analytical approaches [77-79] are also reported. The agreement be- 
tween COM and BA is excellent. The self-consistent Ladder approximation 
(SCLA) [78] shows also a very good agreement for all values of the coupling, 
but it does not have the correct behavior for an infinite value of the Coulomb 
interaction. Moreover, both the Ladder (LA) [78] and the Gutzwiller (GWA 
and GWF) [77] approximations go to zero at finite U, whereas the Renormal- 
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Fig. 1.4. Internal energy E versus: (left) Coulomb interaction JJ for T = and 
n = 1; (right) filling n for T = and = 4. 
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Fig. 1.5. Internal energy £ versus: (left) temperature T for [/ = 4 and n = 0.5 — >■ 1; 
(right) filling n for ?7 = 4 and various values of temperature T. Lanczos data (4 x 4) 
are taken from Ref. [73]. 



ization Group (RG) [79] has the right asymptotic behavior for f/ — ^ oo, but it 
does not reproduce the non-interacting Umit. The doping dependence of the 
internal energy is shown in Fig. 1.4 (right panel) for U = 4. For comparison, 
we also report the BA results [80] and LA and SCLA approaches [78,79]. 
COM agrees reasonably with the BA, reaching the best agreement at half fill- 
ing. The ladder approximation [78] deviates more and more from the BA as 
approaching half filling; the self-consistent ladder approximation [78] probes 
excellently at any doping. In Fig. 1.5 (left panel) we present the internal en- 
ergy versus temperature in the range < T < 5 and [/ = 4 for several values 
of the filling. The results are compared with with finite-temperature Lanczos 
method (FTLM) data [73]. The agreement is quite remarkable in the entire 
range of temperature and for all studied dopant concentrations. In Fig. 1.5 
(right panel), we report the energy per site as a function of the particle density 
for different values of T and U and compare COM results with Lanczos data 
on a 4 X 4 lattice [73] . The agreement is very good on the whole range of filling 
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Fig. 1.6. Specific heat C versus temperature T: (left) at n = 0.5 and U = A; (right) 
at n = 1 and various values of Coulomb interaction U . 

and for all values of temperature. It is worth noticing that at T = 0.1 the nu- 
merical data report a strange kind of oscillations, not present in our results, 
that will be probably absent in the bulk system. It is worth noticing that 
many more comparisons, all showing a very good agreement, between COM 
results and numerical simulation data have been presented in Refs. [81,82]. 
The overall picture emerging from such comparisons is that COM generally 
gives a very accurate description of the behavior of the internal energy as a 
function of the external parameters (T, U and n) over a very wide range of 
their values. 

Specific heat In Fig. 1.6 (left panel), we report the specific heat as a function 
of the temperature for n = 0.5 and U = A. The agreement with FTLM 
data [83] is excellent in the whole range of temperatures. The presence of 
two peaks in the specific heat has been attributed to the spin and charge 
excitations. Our results show a double peak structure too. When U is weak 
the two peaks overlap and there is no resolution. By increasing U the position 
of the charge peak moves to higher temperatures and we distinguish the two 
contributions as shown in Fig. 1.6 (right panel), where C is plotted as a 
function of T at half-filling for several values of U . A peak appears at low 
temperatures when U is rather large {U > 8). This behavior qualitatively 
reproduces the qMC results [84]. It is also relevant to observe the presence 
of a crossing point at a some definite temperature. Such crossing points have 
been experimentally observed in several materials and obtained by means of 
analitical and numerical studies (see Section 3.4.2 of Ref. [5] for references 
and a detailed discussion). 

Entropy In Fig. 1.7 (left panel), the entropy S is reported as a function of 
the filling at [/ = 4 and various values of temperature T. The numerical data 
are taken from Ref. [73]. The agreement between COM and numerical results 
is rather good except for the lower temperature at which quite anomalous 
oscillations appear in the Lanczos data. In Fig. 1.7 (right panel), we report 
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the behavior of the entropy as a function of the temperature in comparison 
with some Lanczos data [73]. The agreement is very good over the whole 
range of temperature and for any reported value of the filling, except for low 
temperatures and low doping. The discrepancy can be related to the capability 
of our data to describe the effects related to the exchange interaction which 
is more and more important as low are the temperature and the doping. 



Two-pole solution - Bosonic Sector 

There is one more way to tackle the problem of computing the response func- 
tions: the two-particle excitations can be considered as a new sector in the 
dynamics of the system and we can choose also for them a suitable basis alike 
it has been done for fcrmions. This approach will be described in this Section. 
We take as bosonic basis the following one [5,42,85] 

= f n.W = ^^.^(^ (1.53) 

This basis is directly generated by the hierarchy of the equations of motion; 
this will assure that the first four bosonic spectral moments have the correct 
functional form [5,8]. The field N^'^^i) obeys the equation of motion 



^Ni^\^) = J(-)W = (-'f.'rt ) (1-54) 



dt ^' ^' \-2dtl^{i) + UK^{i 



where the higher-order bosonic operators are defined as = c^(i)'7p77"(i) — 
r/t(i)a^c°(i)+77t"(i)a^c(i)-ct"(i)a^,7(i) and = c^{)a^r"\i)+c^''\i)<j^c{i)- 
2c^"(i)c7^c"(i) [5,42,85]. We are using the notation c"'(i,i) = Ej"yc(j,t) = 
X^ji ttiiaijcQ, f). The current J^^\i) is projected on the basis (1.53) and the 
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residual current SJ^'^^ is neglected (truncated basis). The normalization and 
energy matrices have the expressions 



4^^(k) ^'^^-u^^^(k)//i^^(k) ; 



4^^(k)//|^^ 

(1.55) 

where ^^^(k) = 4[1 - a(k)]C" and m^^^(k) = -2dt/i^^^(k) + UI,^,^{k). 
The parameter is the electron correlation function C" = (c"(i)ct(i)). The 
quantities /;„p^(k) and /^^p^(k) arc defined as /(„p^(k) = J'([Zp(i, t), p);(j, t)]) 
and /^^^^(k) = J^([K^(i, t), pj^(j, t)]). stands for the Fourier transform op- 
erator. Because of the nonlocality of the bosonic composite field (1.53), the 
analytical form of tliese quantities depends on the dimensionality of the sys- 
tem and it is necessary to separately discuss the different cases [5]. 

The causal G'='('^)(i,i) = {T[N(>''>{i)N^'''>Hj)]} and the correlation C^^^^j) 
{N^^^\i)N^i^^^j)) functions are given by 

G'^^'^Hk.o;) = -2i7rr('')(k)+ (1.56) 



■^a("-^)(k) 



1 + /B(t^) fsjuj) 

w - iot^ (k) +i5 w - iot^ (k) - iS 



(1.57) 



C('^)(k,w) =27rr('^)(k)(5(w) + 27r^5[w-4''Hk)][l + /B(w)](j("'''Hk) . 

(1.58) 



71=1 



The energy spectra a;4'*^(k) are given by wi^'(k) = (-)"-\/e[^^(k)e^^'(k) 
and the spectral functions cr^'^''^\k) have the following expression 

,(7..)(k)=4!W /4^^(k) 1 \ 

The determination of G<='(^)(k,w) and C'^''\k,w) requires the knowledge 
of a set of parameters: correlation functions of fermionic and bosonic oper- 
ators and the unknown F function. According to the self-consistent scheme 
given in Sect. 1.2.6, these parameters are fixed by means of the following 
constraints: (i) all fermionic correlators are calculated within the fermionic 
sector (see previous subsection); (ii) the conservation of the current; (iii) the 
condition that the spin susceptibility be a single value finite function; (iv) the 
local algebra constraint {nij,{i)nij,{i)) = (n) + 2(2^^^o — 1)-D where D is the 
double occupancy. The unknown J'^'^^ (k) function can be calculated either by 
assuming the ergodicity of the system (e.g. /^^^(k) = 6^,fl{2'K/a)'^S^'''>{k)n^) 
or opening a new bosonic sector (the pair sector) [see Ref. [86]]. Once the 
self-consistent equations in the fermionic and in the bosonic sectors have been 
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Fig. 1.8. Spin correlation function S(k) versus momentum k for: (left) [7 = 4, 
T = 0.2 and n = 0.8, 0.33 and 0.19; (right) U = 8,T = 0.2 and T = 0.57 and n = 1, 
0.45 and 0.2. qMC data (8 x 8) are taken from Ref. [87]. 



solved, we can calculate the Green's function and the charge and spin suscepti- 
bilities X/i(l^) = ^ wsll ^ the spin correlation functions 

(1.60) 

Results and comparisons 



correlation function The behavior of 5(k), as a function of the mo- 
mentum, is reported in Fig. 1.8 in comparison with some numerical data [87] 
for different values of filling, Coulomb repulsion and temperature. We have a 
very good agreement with the numerical results for sufficiently high values of 
the doping for all shown values of the Coulomb repulsion. In the proximity of 
half-ffiling the numerical data suffer from a saturation of the antiferromagnetic 
correlation length [76] that becomes comparable with the size of the cluster. 
For [/ = 4 and n = 0.8 (see Fig. 1.8 (left)), the correlation length is slightly 
smaller than the size of the cluster: our solution results capable to describe this 
situation fairly well (the height of the peak at Q is exactly reproduced and the 
momentum dependence is qualitatively correct, again practically exact along 
the diagonal) except for the exact value of the numerical data along the main 
axes. This is probably due again to an overestimation of the correlations by 
the numerical analysis owing to the finite size of the cluster. For U = 8 and 
n = 1 (see Fig. 1.8 (right)), in order to reproduce the numerical data we need 
to increase the temperature as to decrease our value of the correlation length 
and match that of the numerical analysis, which is stuck at the saturation 
value due to the finiteness of the clusters. The results of such a procedure are 
astonishing, we manage to exactly reproduce the numerical data for any value 
of the momentum, and not only at Q, revealing the correctness and power of 
our approach and the limitations of the numerical analysis. 
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Fig. 1.9. Charge correlation function 
A'^(k) versus momentum k for ?7 = 8, 
T = 0.125 and 0.25 and n = 0.155, 0.2 

and 0.5; qMC data (8 x 8, 12 x 12, 16 x 
16) are from Ref. [88]. 
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Fig. 1.10. (left) Spin spectrum ^(k) along the PD for n = 1, T = \QK « O.OOSt 
and U = 8.8f {t k, Q.ieV). (right) The spin-wave intensity along the PD for n = 1, 
T = 295/s: « 0.08t and U = 8.8t [t w 0.3eF). Experimental data are taken from 
Ref. [89]. 



Charge correlation function iV(k) is reported in Fig. 1.9, as a function of the 
momentum, for various fillings and temperatures and ?7 = 8. We have again 
a very good agreement with quantum Monte Carlo results [88] for all shown 
values of the external parameters and of the momentum. The enhancement 
at k = Q = M for n = 0.5 can be interpreted as the manifestation of a 
rather weak ordering of the charge with a checkerboard pattern. COM re- 
sults manage to reproduce such double peak structure showing a capability to 
quantitatively describe, also in a translational invariant phase, rather strong 
charge correlations. 

Comparison with the experimental data for La2CuOi In Fig. 1.10 (left), we 
report the energy spectrum of the spin-spin propagator along the Principal 
Directions and compare with the experimental data of Ref. [89] obtained for 
La2CuOi by means of inelastic magnetic neutron scattering. We have fixed 
the temperature at T = IQK according to the experimental value; the value 
of the transfer integral (t = 0.3eF) and of the Coulomb repulsion {U = 8.8f) 
have been chosen in order to fit the experimental points and they are within 
the ranges {t = 0.3 ± 0.026^; [/ = 2.2 ± QAeV) suggested in Ref. [89]. The 
agreement with the experimental data is very good all over the momentum 
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Fig. 1.11. Static spin susceptibility Xs(k) for various values of temperature T, 
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Fig. 1.12. Static spin susceptibility Xs(k) for [/ = 4, n = 0.97 and T = 0.01 and 
(left) t' = -0.1 [(right) t' = -0.4]. 



space except around (tt, tt). As a matter of fact, the experimental data refer to 
the antiferromagnetic phase of the material [the experimental spectrum gets 
completely soft at (tt, tt)] and our paramagnetic solution obviously cannot fully 
describe such behavior. However, it is worth noting that the Hubbard model 
at half filling presents so strong antiferromagnetic correlations that they also 
show in the paramagnetic phase through a quite deep minimum. The spin- 
wave intensity is reported in Fig. 1.10 (right) and again compared with the 
experimental data. The agreement is again very good over all the momentum 
space and shows once more the capability of the Hubbard model within our 
formulation to catch the physics of such a strongly correlated system. 

Incommensurability In Fig. 1.11, we present the static spin susceptibility 
Xs(k) for various temperatures. The choice for the parameters is U = 4, f = 
—0.1 and the particle filling has been fixed as n = 0.97. t' is the next-nearest- 
hopping integral and the related COM formulation (i.e., for the t-t'-U Hubbard 
model) in the 2-pole approximation can be found in Refs. [40,46-48,95-98]. In 
(left) and (right), x«(k) is reported along the lines k = (tt, fc) and k = {k,k), 
respectively. In both cases by increasing the temperature the incommensurate 
double-peak structure becomes a broad maximum centered at (7r,7r). The 
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Fig. 1.13. (left) Incommensurability amplitude 5{x) versus doping x. COM results 
(solid line) refer to U — 4 and T — 0.01. (right) Experimental values of critical 
temperature Tc{x) for La2-xSrxCu04, taken from Yamada [90], Loram [91], van 
Dover [92] , Shafer [93] and Torrance [94] , versus calculated incommensurability am- 
plitude 5{x). The solid line is a guide to the eyes. 

intensity over the whole Brillouin zone can be clearly observed in Fig. 1.12. 
The important features of the data are: the overall square symmetry of the 
scattering with the sides of the square parallel to the {k, k) and to the (fc, ~k) 
lines and the accumulation of intensity near the corners of the square. These 
features reproduce the experimental situation for La2-x{Ba, Sr)xCu04 [90, 
99 105]. In In Fig. 1.13 (left panel), the incommensurability amplitude S{x) 
is shown as a function of doping. For comparison we report the experimental 
data of Refs. [90, 102, 103, 105]. The linear behavior of ^(a;), observed in the 
low-doping region, agrees exceptionally well with the experimental data. One 
of the most striking features of the results presented in Fig. 1.13 (left panel) 
is the resemblance between the incommensurability amplitude 6{x) and the 
critical temperature Tc. S{x) is maximum in the region of optimal doping 
where is maximum. It has already been experimentally observed in Ref. [90] 
that there is a linear relation between 6{x) and up to the optimal doping 
level X w 0.15. Our theoretical results confirm this behavior and show that a 
close similarity between 6{x) and exists in the entire region of doping. This 
can be seen in Fig. 1.13 (right panel) where experimental values for Tc, taken 
from Refs. [90-94], are reported versus our calculated incommensurability 
amplitude 6{x). 

1.3.2 The residual self-energy i7(k, u>) 

We now come to the problem of taking into account the residual self-energy 
E(k,uj) appearing in the full propagator G'(k, w) once the truncated basis 
(1.30) has been adopted. In particular, we have computed [29-36] X'(k, w) 
in the non-crossing approximation (NCA) by noting that it was possible to 
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Fig. 1.14. Self-consistency sciieme to compute tlie fermionic full propagator G in 
terms of the charge-charge and spin-spin 2-pole propagator B and the residual self- 
energy S in the COM(2p-FNCA(2p)). 



exactly rewrite the residual current SJ{i) as the sums of products of a bosonic 
operator, embodying either charge n{i), spin n(i) or pair p(i) operator, and of 
a fermionic operator, which turned to be either ^{i) or ri{i). Then, following the 
prescription of NCA, we expressed the residual self-energy i7(k, uj) as sums of 
products of the charge-charge and spin-spin propagators (we neglected pair- 
pair propagator) and of the fermonic propagator. The bosonic propagators 
have been computed in the 2-pole approximation as described in the previous 
section. 

The whole framework closes the self-consistency cycle reported in Fig 1.14, 
where Gq and Bq are the fermionic and bosonic (charge and spin) propagators 
in the 2-pole approximation (in parentheses the unknown parameters they de- 
pend on) , 2Jn is the residual self-energy at the n*'* step (it depends on G„ and 
Bn, which are equal to Go and Bq at the startup), G„ is the full fermionic 
propagator (it depends on the unknowns in round parentheses and on Sn-i) 
and, finally, Bn is the bosonic (charge and spin) propagators in the 2-pole 
approximation at the n*^ step (it depends on G„, besides unknown parame- 
ters). The self-consistency cycle terminates when the fermionic parameters /i, 
A, and p do not change; to get six-digit precision for the latter fermionic pa- 
rameters, we usually need about ten cycles (it varies very much with doping, 
temperature, and interaction strength) on a 3D grid of 128 x 128 points in 
momentum space and 4096 Matsubara frequencies. 

In Fig. 1.15 (top panel), the spectral function ^(k,w) along the principal 
directions (r = (0, 0) ^ S* = (7r/2, 7r/2) ^ M = (tt, tt) ^ X = (tt, 0) ^ 5 = 
(7r/2,7r/2) Y — {0,tt) ^ F — (0,0)) is reported for Coulomb repulsion 
U = 8, filling n = 0.92 and temperature T = 0.02 in the frequency range 
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Fig. 1.15. U ^ 8, n = 0.92 and T = 0.02. (top) Spectral function A(k,aj) along 
the principal directions, (bottom left) Spectral function at the chemical potential 
^(k, a; = 0) as a function of momentum k. (bottom right) Imaginary part of the 
self-energy at the chemical potential S"{'k,uj = 0) as a function of temperature 
squared at two momenta: S — (7r/2,7r/2) and 5, where the main diagonal of 
Brillouin zone touches the phantom portion of the hole pocket. 



in the proximity of the chemical potential (uj = 0). We can clearly see may 
unconventional features: kinks in the dispersion along the main diagonal (F — > 
M) and the side of the Brillouin zone (M — > X, Y) at energies comparable with 
the effective exchange energy in the system (J = At'^/U); extended regions 
in momentum where the imaginary part of the self-energy is strong enough 
to selectively suppress the spectral weight; almost doubling of the Brillouin 
zone according to the rather high intensity of the not-so-short-range anti- 
ferromagnetic correlations; formation of a hole pocket centred along the main 
diagonal {F M); formation of almost-closed electron pockets centred at 
X = (tt, 0) and Y — (0, tt); high- intensity of the spectral weight ai X = (tt, 0) 
and Y = (0, tt) (van Hove points) although well below the Fermi surface; no- 
flat dispersion along the anti-diagonal {X Y) generated dynamically {t' is 
not present in the original Hamiltonian); strong suppression of the spectral 
weight in proximity of M = (tTjTt) (again, due to the rather high intensity of 
the not-so-short-range anti-ferromagnetic correlations) which will lead to the 
appearance of a pseudogap. This scenario is just the one claimed for high-Tc 
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cuprates in the underdoped region by ARPES and many other experimental 
techniques. 

In Fig. 1.15 (bottom left panel), the spectral function at the chemical 
potential ^(k, w = 0) as a function of momentum k is reported for Coulomb 
repulsion U = 8, filling n ~ 0.92 and temperature T ~ 0.02. The hole pocket 
and the electron quasi-pocket are now clearly visible. The Fermi surface is 
deconstructed: there is coexistence of a small Fermi surface (the hole pocket) 
and of a large one partially solving the hoary dichotomy signalled by main 
experimental techniques; the actual Fermi surface, the locus formed by the 
momenta where the spectral weight reaches the highest intensity (this also 
the only Fermi surface detectable by ARPES), is open (it is just an arc - half 
of the hole pocket) in contrast to what Fermi liquid picture would require; it 
is just the imaginary part of the residual self-energy to eat up a portion of the 
hole pocket and make this latter just a phantom arc. The residual self-energy 
results stronger and stronger close to M = (7r,7r) as one should expect for 
high-intensity not-so-short-range anti-ferromagnetic correlations. 

In Fig. 1.15 (bottom right panel), the imaginary part of the self-energy 
at the chemical potential S"{k,Lo = 0) as a function of the temperature 
squared at two momenta: S = (7r/2,7r/2) and S. The latter corresponds 
the point where the main diagonal {F — > M) touches the phantom arc of the 
hole pocket; it is clearly visible in Fig. 1.15 both (top panel) and (bottom 
left panel). It is striking the qualitative difference in the dependence of the 
imaginary part of the residual self-energy on in the two points: S lives 
in an ordinary Fermi liquid, where the dependence on temperature is exactly 
quadratic; S belongs to a portion of the system interested to a non-Fermi 
liquid description according to its dependence far from quadratic. Also the 
dependence of frequency (not shown) turn out to be almost perfectly quadratic 
in S and with strong linear components in S_. This unconventional behaviour 
can be at the basis of the anomalous features shown by resistivity in these 
materials. 

1.3.3 Four-pole solution 

In the previous section, we have shown how the 2-pole solution can be im- 
proved by taking into account residual dynamical corrections. However, as 
discussed in Sect. 1.2.5, the 2-pole solution can be also improved by enlarg- 
ing the basis through the addition of higher-order composite operators. This 
aspect is illustrated in this Section. 

As shown in Sect. 1.3.1, the higher order field 7r(i) appears in the equations 
of motion of £,{i) and r]{i). The effect of this field was approximately taken into 
account by projecting it on the basis (2.2). Here we promote this field to the 
rank of new composite operator. Actually, we divide 7r(i) into two operators 
Tr{i) = ^s{i) + Vs{i), similarly to what we have done for the electronic field 
[c{i) = ^{i) + ?7(^)], defined as 
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(1.61) 



and take as new basis 



r]{i) 



(1.62) 



The reason for such a choice is that this basis is the fermionic closed 
basis for the Hubbard model on two sites [9]. In such a case, ^s{i) and r]s{i) 
are cigenoperators of the Hamiltonian with eigenenergies given by E^ih) = 
—jjL + ta{\<i) — Ju and E4^{k) = —ji + ta{\s.) + U + Ju- It is worth noting that 
the energy term Ju = [y/U^~+lQ^ —U]/2 returns, in the limit U/t^l, the 
virtual exchange energy J = At^ /U . 

COM, as formulated in Sect. 1.2, requires as first step the calculation of 
the normalization / and energy e matrices (4x4 matrices). The detailed 
expression of the matrix / is reported in Ref. [17]. In order to effectively 
perform calculations, the elements of / which contain three-site correlation 
functions of the form 



(^"(i)i?"(i)C(i)D(i)) = -{{A{\)B{\)rC{\)D{\)) + 

+ E aijaik(^(j)B(k)C(i)D(i)) (1.63) 



have been computed through the following decoupling procedure: (i) the first 
term is reducible into two-site correlation functions, which can be directly 
evaluated in terms of the propagators under analysis; (ii) the second one 
has been decoupled in terms of two-site correlation functions. This procedure 
preserves the particle-hole symmetry enjoined by the Hamiltonian. 

The energy matrix can be straightforwardly calculated by means of the 
anticommuting algebra given in Tab. 1.1. However, due to the complexity of 
the equations of motion of the fields ^s(z) and r)a{i), a direct calculation will 
lead to the appearance of an enormous number of unknown correlation func- 
tions. In order to determine all correlation functions, we would then be forced 
to use some uncontrolled decoupling procedure and we would completely lose 
every control over the approximation procedure. Hence, we have opted for 
a controlled approximation at the level of equations of motion by neglecting 
irreducible three-site operators and paying attention to evaluating exactly all 
one- and two- site components. By means of this approximation procedure, 
the matrix e can be easily computed. Then, we have all ingredients, following 
the procedure reported in Sect. 1.2.6 and 1.3.1, to calculate all single-particle 
and thermodynamic properties. 
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Fig. 1.16. The dispersion relation at U/t = 8, T/t = 0.5 and (left) n = 0.75 and 
(right) 0.94. The 2-pole solution (dashed hne) and QMC data (circle) of Ref. [106] 
are also reported. 




T/l T/l T/l 



Fig. 1.17. Internal energy E and specific heat C per site at U/t — 4 and 8 and 
(left) n = 0.75 and (right) 0.90. Data from finite temperature Lanczos [73] and 
QMC [84, 107] are provided for comparison. 

Results and comparisons 

In Fig. 1.16, we provide a detailed comparison of the band structure obtained 
by the present formulation with QMC results [106]. As can be easily seen, 
we have a good agreement with the QMC data, especially for the low-energy 
band around the Fermi level. Internal energy E and specific heat C per site at 
U/t = A and 8 and n = 0.75 and 0.90 are reported in Fig. 1.17. Data from finite 
temperature Lanczos [73] and QMC [84, 107] are provided for comparison. As 
regards the internal energy, the agreement with the Lanczos data is excellent 
except for the low temperature region at U/t — 8. As regards the specific 
heat, we observe a sharp peak around T/t = 0.3 and a fairly broad peak 
in the higher temperature region T/t — 1. The two-peak structure is more 
pronounced for U/t ~ 8, but not so much for U/t = 4. This tendency is 
also observed in several numerical simulations [73,84,107,108]. Usually, the 
sharp peak at lower temperatures and the broad peak at higher temperatures 
are interpreted as consequences of spin and charge fluctuations related to 
the energy scales of J and U, respectively. The main difference between our 
results and numerical ones regards the height of the peak in the specific heat 
around T/t = 0.3 that comes from the decrease in the internal energy. This 
is an indication of well established spin ordering which cannot be correctly 
evaluated on a small cluster. Numerical simulation cannot describe spin and 
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charge ordering in the case that the correlation lengths exceed the cluster 
size. On the other hand, in our formulation, there is a tendency to have too 
pronounced spin and charge correlations. 



1.3.4 Superconducting solution 



The various solutions of the Hubbard model presented up to now refer to 
a paramagnetic homogeneous state. However, the same formulation can be 
easily applied to states with completely different symmetry properties. Ferro- 
magnetic, antifcrromagnetic and superconducting solutions of the model have 
been reported in Ref. [5]. As an example, we here summarize the supercon- 
ducting solution of the model in the d-wave channel. 

We consider the Nambu representation and introduce the 4-vector com- 
posite field [5,109,110] 



(1.64) 



At a first level of approximation, the current of this operator is projected 
on the basis and the residual current is neglected. The effects of the residual 
current have also been analysed, but they will not be presented here. Although 
the model admits various superconducting phases with different symmetries 
of the order parameter (see Ref. [5]), on the basis of experimental evidence 
for high superconductors, we restrict the analysis to singlet pairing and 
d-wave symmetry for a two-dimensional lattice. Then, we impose as boundary 
conditions {t{jf{i)tpi{i)) = and {tljf{i)tp'^{i)) = 0. Furthermore, we assume 
translational and spin rotational invariance. Accordingly, the normalization 
and m matrices have the expressions 



7(k) = 



f hi \ 

/22 

111 

\ /22 / 



hi 

I22 



(1.65) 



m = 



nil 

1712 —mi 



mi 



^ / mil mi2 
\^mi2 m22 



m2 = mi3 



-1 1 



(1.66) 



with 



mii(k) = -fihi - 4tA - 4ta(k)(l -n + p) 
TOi2(k) = it A + 4to(k)(p - 122) 
mi3(k) = 4i7(k)/ 

m22(k) = {U- n)l22 - 4:tA - 4ta(k)p 



(1.67) 
(1.68) 
(1.69) 
(1.70) 
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Fig. 1.18. Critical temperature Tc as 
a function of the doping 5 for U = At 
and t ~ G.lTeV . Experimental data for 
La2-sSrsCu04 from Yamada [90], Lo- 
ram [91], van Dover [92], Shafer [93] and 
Torrance [94]. Phenomenological curve 
from Presland [111]. 



and 7(k) = [cos(fcj;a)— cos(fcya)]/2. The parameters A and p have been defined 
in Sect. 1.3.1, the parameter / plays the role of order parameter and is defined 
as / = (c|(i)c^(i)[c;(2)c^(i)]T) + ([c|(i)c^(i)]Tc;(i)c^(i)). The retarded and 
correlation functions, which are now 4x4 matrices containing the anomalous 
components, are then known functions of the internal parameters /.t. A, p 
and /. The self-consistent determination of these parameters is performed by 
means of the equations given in Sect. 1.3.1 and by using a simple decoupling 
procedure for the latter. 

Results and comparisons 

In Fig. 1.18, we report the critical temperature as a function of the doping 
6 for U = At and t « Q.lleV . The experimental data for La2-sSrsCuOi are 
taken from Yamada [90], Loram [91], van Dover [92], Shafer [93] and Tor- 
rance [94]. The phenomenological curve suuggested by Presland [111] is also 
reported. The agreement is really very good: the position and presence of end- 
points (critical fillings at zero temperature), the position of the peak (critical 
filling at highest temperature) and the overall shape. This clearly shows that, 
together with the very remarkable capability to describe the very unconven- 
tional physics of the underdoped region (see Sect. 1.3.2), COM manages to 
catch many relevant aspects of cuprate physics. 

1.4 Conclusions and Outlook 

Developing new methods and techniques to properly take into account electron- 
electron interaction in many-body systems is certainly one of the most chal- 
lenging problems in theoretical physics in the last decades. Among the vari- 
ous methods proposed and reported in this volume, the Composite Operator 
Method is a general method to study strongly correlated systems. The for- 
malism, presented in a systematic way, is based on two main concepts: the 
use of the propagators of composite operators as building blocks of a pertur- 
bation calculation; use of algebra constraints to fix the representation of the 
Green's functions in order to keep the algebraic and symmetry properties of 
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the Hamiltonian. Although the method has been used for the study of several 
models, wc preferred just to present its application to the Hubbard model. 
The presentation has been given by going through a scheme of successive ap- 
proximations, reporting at each stage a comparison with the results of exact 
and numerical data. Some comparison with the experimental data for high 
Tc cuprates has also been reported. It emerges that the Composite Operator 
Method can be considered as a powerful method to obtain information re- 
garding the features of complex strongly correlated Hamiltonian systems and 
materials mimed by them. 
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